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Abstract 

This contribution provides a general framework to use Lagrange multipliers for the simula¬ 
tion of low Reynolds number hber dynamics based on Bead Models (BM). This formalism 
provides an efficient method to account for kinematic constraints. We illustrate, with several 
examples, to which extent the proposed formulation offers a ffexible and versatile framework 
for the quantitative modeling of flexible fibers deformation and rotation in shear flow, the 
dynamics of actuated filaments and the propulsion of active swimmers. Furthermore, a new 
contact model called Gears Model is proposed and successfully tested. It avoids the use of 
numerical artifices such as repulsive forces between adjacent beads, a source of numerical 
difficulties in the temporal integration of previous Bead Models. 

Keywords: Bead Models, hbers dynamics, active hlaments, kinematic constraints, Stokes 
flows 


1. Introduction 

The dynamics of solid-liquid suspensions is a longstanding topic of research while it 
combines difficulties arising from the coupling of multi-body interactions in a viscous fluid 
with possible deformations of flexible objects such as hbers. A vast literature exists on the 
response of suspensions of solid spherical or non-spherical particles due to its ubiquitous 
interest in natural and industrial processes. When the objects have the ability to deform 
many complications arise. The coupling between suspended particles will depend on the 
positions (possibly orientations) but also on the shape of individuals, introducing intricate 
effects of the history of the suspension. 

When the aspect ratio of deformable objects is large, those are generally designated as 
hbers. Many previous investigations of hber dynamics, have focused on the dynamics of rigid 
hbers or rods in El- Compared to the very large number of references related to particle 
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suspensions, lower attention has been paid to the more complicated system of flexible hbers 
in a fluid. 

Suspension of flexible hbers are encountered in the study of polymer dynamics 
whose rheology depends on the formation of networks and the occurrence of entanglement. 
The motion of hbers in a viscous huid has a strong ehect on its bulk viscosity, microstruc¬ 
ture, drainage rate, hltration ability, and hocculation properties. The dynamic response 
of such complex solutions is still an open issue while time-dependent structural changes of 
the dispersed hbers can dramatically modify the overall process (such as operation units in 
wood pulp and paper industry, how molding techniques of composites, water purihcation). 
Biological hbers such as DNA or actin hlaments have also attracted many researches to 
understand the relation between hexibility and physiological properties |5]. 

Flexible hbers do not only passively respond to carrying how gradients but can also be 
dynamically activated. Many of single cell micro-organisms that propel themselves in a huid 
utilize a long hagellum tail connected to the cell body. Spermatozoa (and more generally one- 
armed swimmers) swim by propagating bending waves along their hagellum tail to generate 
a net translation using cyclic non-reciprocal strategy at low Reynolds number [6]. These 
natural swimmers have been modeled by artihcial swimmers (joint microbeads) actuated by 
an oscillating ambient electric or magnetic held which opens breakthrough technologies for 
drug on-demand delivery in the human body [^. 

Many numerical methods have been proposed to tackle elasto-hydrodynamic coupling 
between a huid how and deformable objects, i.e. the balance between viscous drag and 
elastic stresses. Among those, “mesh-oriented” approaches have the ambition of solving a 
complete continuum mechanics description of the huid/solid interaction, even though some 
approximations are mandatory to describe those at the huid/solid interface. Without being 
all-comprehensive, one can cite immerse boundary methods (e.g. IHISIIIDIIII]), extended 
hnite elements (e.g. [I2]), penalty methods [131 E]) particle-mesh Ewald methods [T5|, 
regularized Stokeslets [isiin], Force Coupling Method [T^ . 

In the specihc context of low Reynolds number elastohydrodynamics [IH] , difficulties arise 
when numerically solving the dynamics of rigid objects since the time scale associated with 
elastic waves propagation within the solid can be similar to the viscous dissipation time- 
scale. In the context of self propelled objects the ratio of these time scales is called “Sperm 
number”. When the Sperm number is smaller or equal to one, the object temporal response 
is stih, and requires small time steps to capture fast deformation modes. In this regime, 
huid/structure interaction ehects are difficult to capture. One possible way to circumvent 
such difficulties is to use the knowledge of hydrodynamic interactions of simple objects in 
Stokes how. 

This strategy is the one pursued by the Bead Model (BM) whose aim is to describe 
a complex deformable object by the hexible assembly of simple rigid ones. Such hexible 
assemblies are generally composed of beads (spheres or ellipsoids) interacting by some elastic 
and repulsive forces, as well as with the surrounding huid. For long elongated structures, 
alternative approaches to BM are indeed possible such as slender body approximation m 
|2Tl |22| or Resistive Force Theory [23l [2ll |25] . 

One important advantage of BM which might explain their use among various commu- 
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nities (polymer Physics ESI EH El], micro-swimmer modeling in bio-fluid mechanics 
[Ml EH Sni ESI, flexible flber in chemical engineering [Ml SSI EDI ESj); relies on their para¬ 
metric versatility, their ubiquitous character and their relative easy implementation. We 
provide a deeper, comparative and critical discussion about BM in Section However, we 
would like to stress that the presented model is more clearly oriented toward micro-swimmer 
modeling than polymer dynamics. 

One should also add that BM can be coupled to mesh-oriented approaches in order to pro¬ 
vide accurate description of hydrodynamic interactions among large collection of deformable 
objects at moderate numerical cost [13] • Many authors only consider free drain, i.e no Hy¬ 
drodynamic Interactions (no HI), [271SH SH ED] held interactions associated with the 
Rotne-Prager-Yamakawa tensor [Ml EH EH ES]- This is supported by the fact that far-held 
hydrodynamic interactions already provide accurate predictions for the dynamics of a single 
flexible flber when compared to experimental observations or numerical results. In order to 
illustrate the method we use, for convenience, the Rotne-Prager-Yamakawa tensor to model 
hydrodynamic interactions. We wish to stress here that this is not a limitation of the pre¬ 
sented method, since the presented formulation holds for any mobility problem formulation. 
However, it turns out that for each configuration we tested, our model gave very good com¬ 
parisons with other predictions, including those providing more accurate description of the 
hydrodynamic interactions. 

The paper is organized as follows. First, we give a detailed presentation of the Bead 
Model for the simulation of flexible fibers. In this section, we propose a general formulation of 
kinematic constraints using the framework of Lagrange multipliers. This general formulation 
is used to present a new Bead Model, namely the Gears Model which surpasses existing 
models on numerical aspects. The second part of the paper is devoted to comparisons and 
validations of Bead Models for different configurations of flexible fibers (experiencing a flow 
or actuated filaments). 

Finally, we conclude the paper by summarizing the achievements we obtain with our 
model and open new perspectives to this work. 

2. The Bead Model 

2.1. Detailed Review of previous Bead Models 

The Bead Model (BM) aims at discretizing any flexible object with interacting beads. 
Interactions between beads break down into three categories: hydrodynamic interactions, 
elastic and kinematic constraint forces. Hydrodynamics of the whole object result from 
multibody hydrodynamic interactions between beads. In the context of low Reynolds num¬ 
ber, the relationship between stresses and velocities is linear. Thus, the velocity of the 
assembly depends linearly on the forces and torques applied on each of its elements. Elastic 
forces and torques are prescribed according to classical elasticity theory [55] of flexible mat¬ 
ter. Constraint forces ensure that the beads obey any imposed kinematic constraint, e.g. 
fixed distance between adjacent particles. All of these interactions can be treated separately 
as long as they are addressed in a consistent order. The latter is the cornerstone which 
differentiates previous works in the literature from ours. Numerous strategies have been 
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employed to handle kinematic constraints. 


[321 iol ESI EH and [50] used a linear spring to model the resistance to stretching and 
compression without any constraint on the bead rotational motion (Fig. [^. The resulting 
stretching force reads: 

F" = (1) 

where 


• kg is the spring stiffness, 


• = Fj+i — Fj is the distance vector between two adjacent beads (for simplicity, 
equations and hgures will be presented for beads 1 and 2 and can easily be generalized 
to beads i and i + 1), 

• r? 2 vector corresponding to equilibrium. 


However, regarding the connectivity constraint, the spring model is somehow approxi¬ 
mate. A linear spring is prone to uncontrolled oscillations and the problem may become 
unstable. Many other authors, among which [2HI EHl EO], thus use non-linear spring models 
for a better description of polymer physics. Nevertheless, the repulsive force stiffness has an 
important numerical cost in time-stepping as will be discussed in section 2.6.3 Furthermore, 
unconstrained bead rotational motion leads to spurious hydrodynamic interactions and thus 
limits the range of applications for these BM. 


Alternatively, [191 mUMl HZ] and US] constrained the motion of the beads such that 
the contact point for each pair q remains the same. While more representative of a flexible 
object, this approach exhibits two main drawbacks: 

1. a gap between beads is necessary to allow the object to bend (see Fig. [^, 

2. it requires an additional center to center repulsive force, and thus more tuning numer¬ 
ical parameters to prevent overlapping between adjacent beads. 

Consider two adjacent beads, with radius a, linked by a hinge Ci (typically called ball 
and socket joint). The gap Eg dehnes the distance between the sphere surfaces and the joint 
(see Fig. [^. Denote p* the vector attached to bead i pointing towards the next joint, i.e. 
the contact point q. 

The connectivity between two contiguous bodies writes: 

[ri (a -F £g)pi] - [r2 - (a -h e:g)p2] = 0 ( 2 ) 

and its time derivative 

[fi - (a + eg)pi X u)i] - [f2 + {a + eg)p2 x 0^2] = 0. ( 3 ) 

f j and (jJi are the translational and rotational velocities of bead i. The constraint forces and 
torques associated to ([^ are obtained either by solving a linear system of equations involv¬ 
ing beads velocities [S3]) or by inserting (i into the equations of motion when neglecting 
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Figure 1: Spring model: linear spring to keep constant the inter-particle distance. 



Figure 2: Joint Model: overlapping due to bending if no gap between beads. 



Figure 3: Joint Model: Ci is separated by a gap Sg from the beads. 



Figure 4: Gears Model: contact velocity must be the same for each bead (no-slip condition). 
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hydrodynamic interactions [IHl EH] . 

The gap width 2eg controls the maximnm cnrvature allowed withont overlapping. From 
the sine rnle, one can derived the simple eqnation relating Eg and 



Once aware of these limitations, the gap Eg, range and strength of the repnlsive force shonld 
be prescribed depending on the problem to be addressed. 

[56] and |13| proposed a more sophisticated Joint Model than those hitherto cited, nsing a 
fnll description of the links dynamics along the cnrvilinear abscissa. They derived a snbtle 
constraint formnlation which ensnres that the tangent vector to the centerline is continnous 
and that the length of links remains constant. These two works are worth mentioning since 
they avoid an empirical tnning of repnlsive forces. Yet, [56] compnted the constraint forces 
and torques with an iterative penalty scheme instead of using an explicit formulation. 


Finally, it is worth mentioning that the bead model proposed in [3T] circumvents the 
inextensibility difficulty by imposing constraints on the relative velocities of each successive 
segments, so that their relative distance is kept constant. Using bending potential, [21] 
permit overlap between beads with restoring torque (cf. Fig. |^. A Lagrangian multiplier 
formulation of tensile forces is also used in isii, which is equivalent to a prescribed equal 
distance between successive beads. Again, inextensibility condition does not prevent bead 
overlapping due to bending in this formulation. The computation of contact forces which 
is proposed in the following section 2.2 generalizes the Lagrangian multiplier formulation 
of [21] to generalized forces. Using more complex constraints involving both translational 
and angular velocities, we show that it is possible to accommodate both non-overlapping 
and inextensibility conditions without additional repulsive forces (using the rolling no-slip 
contact with the gears model detailed in 2.3). This proposed general formulation is also well 


suited for any type of kinematic boundary conditions as illustrated in Section 3.4 


2.2. Generalized forces, virtual work principle and Lagrange multipliers 

The model and formalism proposed in this article rely on earlier work in Analytical Me¬ 
chanics and Robotics [HE] |59] . The concept of generalized coordinates and constraints which 
has proven to be very useful in these contexts is described here. Generalized coordinates 
refer to a set of parameters which uniquely describes the conhguration of the system relative 
to some reference parameters (positions, angles,...). For describing objects of complex shape, 
let us consider the position rj of each bead i G {1, with associated orientation vector 
which is dehned by three Euler angles p = {9, 0, fj). In the following, any collection of vector 
population (ri, ..r*, ..r^vj = R will be capitalized, so that R is a vector in Hence the 

collection of orientation vectors pi will be denoted P, which is a vector of length 3Nb, the 
collection of velocities ^ = Vj, will be denoted V, the collection of angular velocity 

Pi = Ui will be f2, the collection of forces fj, F, the collection of torques 7 ^, F. All V, fl, F 
and F are vectors in 
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Let us then define some generalized coordinate for each bead, which is defined by q* = 
(rj, Pi) = {ri j, r2,i, r^^i, 9i, 'ijji} so that the collection of generalized positions (qi, ..q*, ..qArJ = 
Q is a vector in Generalized velocities are then defined by vectors q* = (vj,LJj) with 

associated generalized collection of velocities Q. 

Articulated systems are generically submitted to constraints which are either holonomic, 
non-holonomic or both [SS]- Holonomic constraints do not depend on any kinematic param¬ 
eter (i.e any translational or angular velocity) whereas non-holonomic constraints do. 

In the following we consider non-holonomic linear kinematic constraints associated with 
generalized velocities of the form [60] 


^Q + B = 0, (5) 

such that ^ is a 6Ni, x A), matrix and B is a vector of A), components. N^. is the number of 
constraints acting on the beads. B and ^ might depend (even non-linearly) on time t 
and generalized positions Q, but do not depend on any velocity of vector Q, so that relation 
(j^ is linear in Q. In subsequent sections, we provide specific examples for which this class 
of constraints are useful. Here we describe, following [UUl 135] how such constraints can be 
handled thanks to some generalized force that can be defined from Lagrange multipliers. The 
idea formulated to include constraints in the dynamics of articulated systems is to search 
additional forces which could permit to satisfy these constraints. First, one must rely on 
generalized forces fi = (fi,7j) which include forces and torques acting on each bead, whose 
collection (fi, fi, ..fArJ is denoted Generalized forces are defined such that the total work 
variation 5W is the scalar product between them and the generalized coordinates variations 
5Q 

6W = d ■SQ = F -SR + r ■ 6P, ( 6 ) 

so that, on the right hand side of ([^ one also gets the translational and the rotational 
components of the work. Then, the idea of virtual work principle is to search some virtual 
displacement 5Q that will generate no work, so that 

d-5Q = 0. ( 7 ) 

At the same time, by rewriting ([^ in differential form 

^dQ -|- Rdt = 0 , ( 8 ) 

admissible virtual displacements, i.e those satisfying constraints ([^, should satisfy 

^5Q = 0. (9) 

Gombining the N^, constraints ([^ with ([^ is possible using any linear combination of these 
constraints. Such linear combination involves Nc parameters, the so-called Lagrange multi¬ 
pliers which are the components of a vector A in Then from the difference between ([^ 
and the Nc linear combination of ([^ one gets 


(d - X- ^)-dQ = 0. 
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( 10 ) 


Prescribing an adequate constraint force 

5c = A-^, (11) 

permits to satisfy the required equality for any virtual displacement. Hence, the con¬ 
straints can be handled by forcing the dynamics with additional forces, the amplitude of 
which are given by Lagrange multipliers, yet to be found. Note also, that this first result im¬ 
plies that both translational forces and rotational torques associated with the Nc constraints 
are both associated with the same Lagrange multipliers. 

This formalism is particularly suitable for low Reynolds number flows for which trans¬ 
lational and angular velocities are linearly related to forces and torques acting on beads by 
the mobility matrix M 



yoo 


+ C : E“. 


( 12 ) 


= (v“,...,v“J and fl°° = ..., correspond to the ambient flow evaluated at 

the centers of mass r*. is the rate of strain 3x3 tensor of the ambient flow. C is a 
third rank tensor called the shear disturbance tensor, it relates the particles velocities and 
rotations to E“ |5l]. Matrix M (and tensor C) can also be re-organized into a generalized 
mobility matrix ^ (generalized tensor ^ resp.) in order to define the linear relation between 
the previously defined generalized velocity and generalized force 


Q = ^5 + ^ ; E° 


( 13 ) 


where ..., . The explicit correspondence between the classical ma¬ 

trix M and the hereby proposed generalized coordinate formulation ^ is given in Appendix 


[A] Hence, as opposed to the Euler-Lagrange formalism of classical mechanics, the dynam¬ 
ics of low Reynolds number flows does not involve any inertial contribution, and provide a 
simple linear relationship between forces and motion. In this framework, it is then easy to 
handle constraints with generalized forces, because the total force will be the sum of the 
known hydrodynamic forces '^h, elastic forces '^e, inner forces associated to active hbers '^a 
and the hereby discussed and yet unknown contact forces to verify kinematic constraints 


d = + with (14) 

+ + (15) 


Hence, if one is able to compute the Lagrange multipliers A, the contact forces will provide 
the total force by linear superposition (14), which gives the generalized velocities with (13). 
Now, let us show how to compute the Lagrange multiplier vector. Since the generalized 
force is decomposed into known forces and unknown contact forces = A ■ ^, relations 
(14) and (13) can be pooled together yielding 

= = Qi- -.YT. (16) 
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So that, using ([^, 

= -B - ^ + ^ : E“), (17) 

one gets a simple linear system to solve for finding A, where stands for the transposition 
of matrix ^. 

2.3. The Gears Model 

The Euler-Lagrange formalism can be readily applied to any type of non-holonomic con¬ 
straint such as (j^. In the following, we propose an alternative model based on no-slip 
condition between the beads: the Gears Model. This constraint, hrst introduced in a Bead 
Model (BM) by [2Zj, conveniently avoid numerical tricks such as artihcial gaps and repulsive 
forces. 

However, 1271 and [61] relied on to an iterative procedure to meet requirements. Here, we use 
the Euler-Lagrange formalism to handle the kinematic constraints associated to the Gears 
Model. 

Gonsidering two adjacent beads (Fig. |^, the velocity Vc^ at the contact point must be 
the same for each sphere: 


Cl 


< = 0 . 


(18) 


and are respectively the rigid body velocity at the contact point on bead 1 and bead 


2. Denote cr^ the vectorial no-slip constraint. (18) becomes 


i.e. 


(ri,a;i.r2,uJ2) ^ 0, 


[r, - aei 2 X uj,] - [ 1-2 - 0621 X UI2] — 0, 


(19) 


( 20 ) 


where ei 2 is the unit vector connecting the center of bead 1 , located at ri, to the center 
of bead 2, located at r 2 (ei 2 = 62 — 61 ). The orientation vector attached to bead i, is 
not necessary to describe the system. Hence, from (20) one realises that cr^ is linear in 


translational and rotational velocities. Therefore equation (19) can be reformulated as 


(tMQ) =j‘Q = 0. 


( 21 ) 


where, Q is the collection vector of generalized velocities of the two-bead assembly 

Q = [fi, CJI, 1-2, n;2]^, (22) 

is the Jacobian matrix of cr^: 

dal 


4 = 5^, 4 = 1,...,3,1 = 1,...,12, 


(23) 
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(24) 


l3 aej, ] , 


and 



(25) 


For an assembly of Nf, beads, iV^ — 1 no-slip vectorial constraints must be satisfied. The 
Gears Model (GM) total Jacobian matrix is block bi-diagonal and reads 


f J 2 


^gm ^ 




J2 J2 


V 


(26) 


-rVi,-! -wNi, — l 


where is the 3x6 Jacobian matrix of the vectorial constraint a for the bead (3. 
The kinematic constraints for the whole assembly then read 


_/0"Q = 0, 


(27) 


The associated generalized forces '^c are obtained following Section 2.2 


2.4- Elastic forces and torques 

We are considering elastohydrodynamics of homogeneous flexible and inextensible hbers. 
These objects experience bending torques and elastic forces to recover their equilibrium 
shape. Bending moments derivation and discretization are provided. Then, the role of 
bending moments and constraint forces is addressed in the force and torque balance for the 
assembly. 

2.4- 1- Bending moments 

The bending moment of an elastic beam is provided by the constitutive law 


m(s) = K^t X 

ds 


(28) 


where K^{s) is the bending rigidity, t is the tangent vector along the beam centerline 
and s is the curvilinear abscissa. Using the Frenet-Serret formula 


dt 

ds 


= Kn, 


(29) 


m(s) = K'^nh, 
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the bending moment writes 


(30) 





75 = -m(s4) 



Figure 5: Beam discretization and bending torques computation of beads 1, 3 and 5. Re¬ 
maining torques are accordingly obtained: 'y\ = m (ss) and 'y\ = —m (ss). 


where k{s) defines the local curvature, n(s) and b(s) are the normal and binormal vectors of 
the Frenet-Serret frame. When the link considered is not straight at rest, with an equilibrium 
curvature (30) is modihed into 

m(s) = {k — b. (31) 


Here, the beam is discretized into W — 1 rigid rods of length I = 2a (Cf Fig. |^. 
Inextensible rods are made up of two bond beads and linked together by a flexible joint 
with bending rigidity K^. Bending moments are evaluated at joint locations Si = {i — 1)1 
for i = 2,N}, — 1, where Sj correspond to the curvilinear abscissa of the mass center of 
the bead. 

The bending torque on bead i is then given by 


7 - = m (si+i) - m (si_i), (32) 

with m (sj) = K^k (sj) b (sj). See Fig. [^for the torque computation on a beam discretized 
with four rods. 


The local curvature k (sj) is approximated using the sine rule |12] 


K(Si) 


1 /T+”e~)~e~^ 

aV 2 


(33) 


where is the unit vector connecting the center of mass of bead i —1 to the center of mass 
of bead i. This elementary geometric law provides the radius of curvature R{si) = 1/^(5*) 
of the circle circumscribing neighbouring bead centers rj_i, r* and rj+i. 

A more general version of the discrete curvature proposed in [HS] can also be used in the 
case of three dimensional motion. In that case, the curvature of the hber is discretized as 
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in 


/ \ ^ 2 — 1,2 ^ ^ 2 , 2+1 

=- 2 ^- 

where, again, ej_i_j is the unit vector connecting the center of mass of bead i — 1 to the 
center of mass of bead i. The bending moment reads 

m (sj) = (si). ( 35 ) 



To include the effect of torsional twisting about the axis of the hber, one would have to 
compute the relative orientation between the frames of reference attached to the beads 
using Euler angles [SB] (see Section 2.2) or unit quaternions as in [S3]. This would provide 


the rate of change of the twist angle along the hber centerline and thus the twisting torque 
acting on each bead. In the following, only bending effects are considered. 


2.4-2. Force and moment acting on each head 

The Gears Model proposed in this paper does not need to consider gaps to allow bending, 
also ensures the connectivity condition and circumvent the use of repulsive forces as 
distances between adjacent bead surfaces remain constant. More specihcally, the tangential 
components of the force Fc, which is only one part of the generalized force 'Sc-, acts as tensile 
force. 

For each bead i, contact forces applied from bead i to bead i + 1 at contact point q 
between bead i and z + 1 (Figure]^ for two beads) is denoted From Newton third law at 
contact point Cj, the contact force applied to bead i from bead z + 1 is obviously —fc^- Total 
force acting on bead z from contact, and hydrodynamic forces f/* reads 


f* = fc,_i - fc, + (36) 

Similarly, the contact force at point Cj produces a moment rric. = atj x f^. associated with 
local tangent vector tj = and distance a to the neutral hber at point c*. Total moment 
acting on bead z from contact points moments, elastic and hydrodynamic torques are then 


li = m, 


Gi —1 


— 


I't+ii 


(37) 


The contribution of contact force and contact moment acting on bead z exactly equals 
the contribution of the generalized contact force. Indeed, using the kinematic constraints 
Jacobian (26) in ([II|, and computing the force and torque contributions, one exactly recovers 


the hrst and the second contributions of the right-hand-side of (36) and (37). In Appendix 


[B] we also show that this model is consistent with classical formulation for slender body 
force and moment balance when the bead radius tends to zero. 


2.5. Hydrodynamic coupling 

Moving objects (rigid or hexible hbers) in a viscous huid experience hydrodynamic forc¬ 
ing. The interactions are mediated by the huid how perturbations which can alter the motion 
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and the deformation of the hbers in a moderately concentrated suspension. The existence 
of hydrodynamic interactions has also an effect on a single hber dynamics while different 
parts of the hber can respond to the ambient how but also to local how perturbations re¬ 
lated to the hber deformation. Resistive Force Theory (RFT) can be used to estimate the 
hber response to a given how assuming that the hber is modeled by a large series of slender 
objects [nU 123] • Slender body theory has also been used [201 ES] to relate local balance 
of drag forces with the hlament forces upon the huid resulting in a dynamical system to 
model the deformation of the hber centerline. This model provided interesting results on 
the stretch-coil transition of hbers in vortical hows. 

In our beads model, the hber is composed of spherical particles to account for the hnite 
width of its cross-section. The hydrodynamic interactions are provided through the solution 
of the mobility problem which relates forces, torques to the translational and rotational ve¬ 
locities of the beads. This many-body problem is non-linear in the instantaneous positions 
of all particles of the system. Approximate solutions of this complex mathematical problem 
can be achieved by limiting the mobility matrices to their leading order. The simplest model 
is called free drain as the mobility matrix is assumed to be diagonal neglecting the HI with 
neighbouring spheres. Pairwise interactions are required to account for anisotropic drag 
effects within the beads composing the hber. The Rotne-Prager-Yamakawa (RPY) approxi¬ 
mation is one of the most commonly used methods of including hydrodynamic interactions. 
This widely used approach has been recently updated by Wajnryb et ah [53] for the RPY 
translational and rotational degrees of freedom, as well as for the shear disturbance tensor 
C which gives the response of the particles to an external shear how (12). 


2.6. Numerical implementation 
2.6.1. Integration scheme and algorithm 

The kinematics of the constrained system results from the superposition of individual 
bead motions. Positions are obtained from the temporal integration of the equation of 
motion with a third order Adams-Bashforth scheme 

dYi 



where r*, Vj are the position and translational velocity of bead i. 

The time step At used to integrate (|3^ is hxed by the characteristic bending time [16| 


At < 


/i(2a)^ 


(39) 


where /i is the suspending huid viscosity. 

The evaluation of bead interactions must follow a specihc order. Elastic and active forces 
can be computed in any order. Constraint forces and torques must be estimated afterwards 
as they depend on 'S'. Then velocities and rotations are obtained from the mobility relation. 
And hnally, bead positions are updated. 
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• Initialization: positions rj(0), 


• Time Loop 


1 . 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 


Evaluate mobility matrix ^ (Q) and (see Section 2.5), 

Calculate local curvatures (33) and bending torques 7 ^ (32) to get '^e, 

Add active forcing and ambient velocity if any, 

Compute the Jacobian matrix associated with non-holonomic constraints ^ (Q), 
Solve ( [TT] ) to get the constraint forces dc = , 

Sum all the forcing terms d = de + da + dc, 

Apply mobility relation (13) to obtain the bead velocities Q, 

Integrate (38) to get the new bead positions. 


2.6.2. Implementation of the Joint Model 

To provide a comprehensive comparison with previous works, we exploit the flexibility 
of the Euler-Lagrange formalism to implement the Joint Model as described in [IH] supple¬ 
mented with hydrodynamic interactions. The joint constraint for two neighbouring beads 
reads 


Itj - (a + iTjjpi X u>j| - |rj+i + (a + <r,)Pi+i x u)j+,] = 0 , 


(40) 


Using the Euler-Lagrange formalism, (40) is reformulated with the Joint Model (JM) 
Jacobian matrix 


= 0 , 


where has the same structure as in (26) and 


(41) 


j- = [j; j‘] 

= [ I3 -(o + ir5)p'' -I3 -(o + £5)p?+i ] . 


(42) 


Accordingly, the corresponding set of forces and torques dc are obtained from Section 


when bending occurs. A repulsive force is added according to |1B] (the force prohle 
proposed by |19] is very stiff, thus very constraining for the time step): 


2.2 As mentioned in Section such formulation does not prevent beads from overlapping 


FL 


' 

-Fo exp 



f _dij + JdA 


dij 



■"O’ 


dij < —Sd, 


—6d < dij < 6 d, 
Tij > 6 d- 


(43) 


df) is an artihcial surface roughness, dij is the surface to surface distance, dij < 0 
indicates overlapping between beads i and j. do is a numerical damping distance which has 
to be tuned to prevent overlapping. Fq is the repulsive force scale chosen in order to avoid 
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numerical instabilities. To deal with this issue, [10] proposed to evaluate Fq from bending 
and viscous stresses. A slight modihcation of their formula for inertialess particles yields 


Fo = CiQu^L (v°° — v) + (^2^ 


I RbE’^ 


(44) 


the bar denotes the average over the constitutive beads or joints where Ci and C 2 are 
adjustable constants. is the bending energy 


Nb-l 

i=l 


(45) 


Bending moments are evaluated at the joint locations s/ = {a+eg) + {i — l) x2{a+eg), i = 
1,Nb — 1. Joint curvature is given by 


n(s-) = 


a + e. 


1 + p*.pi+i 


Similarly to (32), bending torque on bead i is 


= m {si) - m (s/_i) . 

Bead orientation pj is integrated with a third order Adams-Bashforth scheme 

dpi 


dt 


= UJiX Pi. 


(46) 


(47) 


(48) 


The procedure is similar to the Gears Model, pj are initialized together with the posi¬ 
tions. The repulsive force is added to 'S' and can be computed between step 1 and 5 of 


the aforementioned algorithm. Time integration of (48) is performed at step 8 


2.6.3. Constraints and numerical stability 

At each time step, the error on kinematic constraints e is evaluated, after application of 


the mobility relation (13), between step 7 and step 8: 


1/2 


2=1 


for the Gears Model, and 

e^^(t) = 

for the Joint Model. 




(49) 


(60) 


To verify the robustness of both models and Lagrange formulation, a numerical study is 
carried out on a stiff conhguration. 
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A fiber of aspect ratio Vp = 10 with bending ratio BR = 0.01 is initially aligned with a shear 
flow of magnitude 7 = 5s“^. For this aspect ratio, Nf, = 10 beads are used to model the 
hber with the Gears Model. 


Joint Model involves additional items to be hxed. Nb = 9 spheres are separated by a gap 
width 2eg = 0.25a. The repulsive force is activated when the surface to surface distance dij 
reaches the artihcial surface roughness Sjy = 2{a + eg)/19. The remaining coefficients are set 
to reduce numerical instabilities without affecting the Physics of the system: do = (a + £g)/4, 
Cl = 5 and C 2 = 0.5. 

Fig. @ shows the evolution of the maximal mean deviation from the no-slip/joint con¬ 
straint Im = max e{t)/{Nb — 1) normalized with the maximal shear velocity jL depending 

on the dimensionless time step 7 At. First, one can observe that for both Joint and Gears 
models, ImUL weakly depends on 7 At and the resulting motion of the beads complies very 
precisely with the set of constraints, within a tolerance close to unit roundoff (< 2 . 10 “^®). 
Secondly, Joint Model is unstable for time steps 100 times smaller than Gears Model. The 
onset for numerical instability indicates that the repulsive force stiffness dominates over 
bending, thus dictating and restricting the time step. 

As a comparison, [16] matched connectivity constraints within 1 % error for each hber seg¬ 
ment. To do so, they had to use an iterative scheme reducing the time step by 1/3 each 
iteration to meet requirements and limit overlapping between adjacent segments. For similar 
results, a stiff conhguration, such as the sheared hber, is therefore more efficiently simulated 
with the Gears Model. 

Thirdly, inset of Fig. |^shows that, for a given time step, the Gears Model constraints Im/iL 
are satished whatever the shear magnitude. Hence, (39) ensures unconditionally numerical 


stability as bending is the only limiting ehect for the Gears Model. 

Hence, the robustness of the Euler-Lagrange formalism and the numerical integration we 
chose provide a strong support to the Gears Model over the Joint Model. 


As a hnal remark to this section, it is important to mention that the numerical cost 
of the proposed method strongly depends on the choice for the mobility matrix computa¬ 
tion, as usual for bead models. If the mobility matrix is computed taking into account 
full hydrodynamic interactions with Stokesian Dynamics, most of the numerical cost will 
come from its evaluation in this case. This limitation could be overcame using more sophis¬ 
ticated methods such as Accelerated Stokesian Dynamics [HS] or Force Goupling Method 
|18j . Moreover, when considering Rotne-Prager-Yamakawa mobility matrix, its cost only 
requires the evaluation of O {{QNbY) terms. Furthermore, the main algorithmic complexity 
of bead models does not come from the time integration of the bead positions which only 


requires a matrix-vector multiplication (13) at a O {{QNbY) cost. Fast-multipole formulation 
of a Rotne-Prager-Yamakawa matrix can even provide a O {QNb) cost for such matrix-vector 
multiplication |67j . 

The main numerical cost indeed comes from the inversion of the contact forces problem 
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Figure 6 : Dependence of the constraints eu/lL on the time step + : Gears Model, 


□ : Joint Model. Inset: with the Gears Model for a hxed time step given by (39) 

for different values of 7 . 


(HD- It is worth noting that this linear problem is NcXNc which is slightly different from N^x 
Nh, but of the same order. Furthermore, problem ( [ItI ) gives a direct, single step procedure 
to compute the contact forces, as opposed to previous other attemps iZl SSI EEl which 
required iterative procedures to meet forces requirements, involving the mobility matrix 
inversion at each iteration. The cost for the inversion of (17) lies in-between 0{N^) and 
0{N^) depending on the inversion method. 


3. Validations 

3.1. Jeffery Orbits of rigid fibers 

Much of our current understanding of the behavior of hbers experiencing a shear flow has 
come from the work of Jeffery |SH] who derived the equation for the motion of an ellipsoidal 
particle in Stokes flow. The same equation can be used for the motion of an axisymmetric 
particle by using an equivalent ellipsoidal aspect ratio. Rigid hbers can be approximated by 
elongated prolate ellipsoids. An isolated hber in simple shear how rotates in a periodic orbit 
while the center of mass simply translates in the how (no migration across streamlines). The 


period T (51) is a function of the aspect ratio of the hber and the how shear rate while the 


orbit depends on the initial orientation of the object relative to the shear plane 

27r(re -h 1/re) 


T = 


7 


(51) 


17 















100 



Figure 7: Tumbling period T depending on fiber aspect ratio Vp. 


- : theoretical law (51) 

given by (53), —B- : theoretical law (51) with given by (52), A 


with r, 

V : Joint Model. 


Gears Model, 


7 is the shear rate of the carrying flow. r^. is the equivalent ellipsoidal aspect ratio which 
is related to the fiber aspect ratio (length of the fiber over diameter of the cross-section 
which turns out to Vp = Nb with Nb beads). The fiber is initially placed in the plane of 
shear and is composed on Nb beads. No gaps between beads is required in the Joint Model 
because the fiber is rigid and flexibility deformations are negligible. We have compared the 
results with two relations for Ve'- Cox |1] 


1.24rp 

and Larson [69] 


(52) 


r, = 0.7r,. (63) 

This classic and simple test case has been successfully validated in [271 ESI El] • Both the 
Joint and Gears models give a correct prediction of the period of Jeffery orbits (Fig. [^. The 
scaled period yT of simulations remains within the two evolutions based on equations 


and (53). We have tried to compare with the linear spring model proposed by Gauger and 
Stark [10] (and used by Slowicka et ah [00] with a more detailed formulation of hydrodynamic 
interactions). In this latter model, there is no constraint on the rotation of beads and the 
simulations failed to reproduce Jeffery orbits (the fiber does not flip over the axis parallel 
to the flow). 
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3.2. Flexible fiber in a shear flow 

The motion of flexible fibers in a shear flow is essential in paper making or composite 
processing. Prediction and control of fiber orientations and positions are of particular in¬ 
terest in the study of flocks disintegration. Many models have been designed to predict 
fiber dynamics and much experimental work has been conducted. The wide variety of fiber 
behaviors depends on the ratio of bending stresses over shear stress, which is quantified by 
a dimensionless number, the bending ratio BR [701IS3] 


BR 


£'(ln 2re — 1.5) 


(54) 


E is the material Young’s Modulus and /x is the suspending fluid viscosity. 

In the following, we investigate the response of the Gears Model with known results on 
flexible hber dynamics. 


3.2.1. Effect of permanent deformation 

noiin] analysed the motion of flexible threadlike particles in a shear flow depending on 
BR. They observed important drifts from the Jeffery orbits and classihed them into cate¬ 
gories. Yet, the goal of this section is not to carry out an in-depth study on these phenomena. 
Instead, the objective is to show that our model can reproduce basic features characteristic 
of sheared flexible filaments. We analyze first the influence of intrinsic deformation on the 
motion. 

If a fiber is straight at rest, it will symmetrically deform in a shear flow. When aligned with 
the compressive axes of the ambient rate of strain E°°, the fiber adopts the “S-shape” ob¬ 
served in Fig. 8a When aligned with the extensional axes, tensile forces turn the rod back 
to its equilibrium shape. This symmetry is broken when the hlament is initially slightly 
deformed or has a permanent deformation at rest, i.e. a nonzero equilibrium curvature 

> 0 . An initial small perturbation of the shape of a straight hlament aligned with how 
can induce large deformations during the orbit. This phenomenon is known as the buckling 
instability whose onset and growth are quantihed with BR [Z21E3|. Fig. illustrates the 


evolution of a hexible sheared hlament with BR = 0.04 and a very small intrinsic deforma¬ 
tion = 1/(100L). The equilibrium dimensionless radius of curvature is = 200. 

During the tumbling motion it decreases to a minimal value of 2Rmin/L = 0.26. Buckling 
thus increases by 770 times the maximal hber curvature. 


3.2.2. Maximal fiber curvature 

m measured the radius of curvature R of sheared hber for aspect ratios Vp ranging from 
283 to 680. They reported on the evolution of the minimal value Rmini he. the maximal 
curvature Umax, with BR. [53] used the Joint Model with prolate spheroids but no hydro- 
dynamic interactions and compared their results with m- Both experimental results from 
m and simulations from [53] are accurately reproduced by the Gears Model. 


Hydrodynamic interactions between hber elements play an important role in the bending 
of hexible hlaments [731 SSI ED]. As mentioned in Section |2.5| the use of spheres to build 
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xj L 

(a) (b) 

Figure 8 : Orbit of a flexible filament in a shear flow with BR = 0.04. Temporal evolution 
is shown in the plane of shear flow, a) Symmetric “S-shape” of a straight filament, = 0. 
b) Buckling of a permanently deformed rod with an intrinsic curvature = 1/(100L). 


any arbitrary object is well suited to compute these hydrodynamic interactions. However, 
modeling rigid slender bodies in a strong shear flow becomes costly when increasing the fiber 
aspect ratio. First, the aspect ratio of a fiber made up of Ni, spheres is Vp = Ni,. Each time 
iteration requires the computation of ^ and ^ and the inversion of a linear system 

corresponding to Nc relations of constraints with Nc> 3 {Nh — 1 ). Secondly, for a given 


shear rate 7 and bending ratio BR, the Young’s modulus increases as r^. According to p9|). 


the time step becomes very small for large E. [53] partially avoided this issue by neglecting 
pairwise hydrodynamic interactions is diagonal), and by assembling prolate spheroids 
of aspect ratios re ~ 10 . 


Yet, it is shown in Fig. 9a that for a fixed BR, 2Rmin/L converges asymptotically to a 


constant value with Vp. An asymptotic regime (relative variation less than 2 %) is reached 
for Tp > 25. Choosing Vp = 35 thus enables a valid comparison with previous results. 


Our simulation results compare well with the literature data (Fig. 9b) and better match 
with to experiments than [53] . When BR > 0.04, the Gears Model clearly underestimates 
measurements for = 1/(10L) and overestimates them for = 0 . However, Salinas 
and Pittman m indicated that the error quantification on parameters and measurements is 
difficult to estimate as the fibers were hand-drawn. Notably, drawing accuracy decreases for 
large radii of curvature, which leads to the conclusion that the hereby observed discrepancy 
might not be critical. They did not report the value of permanent deformation for the 
fibers they designed, whereas, as evidenced by CB, it has a strong impact on Rmin- A 
numerical study of this dependence should be conducted to compare with m!, Fig. 7. 

used the same approach as [Ml with hydrodynamic interactions to repeat numerically 
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Figure 9: a) Minimal radius of curvature depending on fiber length for several bending 
ratios. ^ : BR = 0.01, ^ : BR = 0.03, ; BR = 0.04 , : BR = 0.07. b) 

Minimal radius of curvature along BR. o ; current simulations with aspect ratio Vp = 35 
and intrinsic curvature = 0 ; • : current simulations with aspect ratio Vp = 35 and 
intrinsic curvature = 1/(10L) ; simulation results from [53] with = 1/(10L): ( 0 : 
Tp = 50, A : Vp = 100, V '■ Tp = 150, □ : rp = 280) ; + : experimental measurements 
from [71|, Tp = 283. 
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the experiments from m ; but their results, though reliable, were displayed such that direct 
comparison with previous work is not possible. 

To conclude, it should be noted that, in |53], the aspect ratio does not affect 2Rmin/L for a 


fixed BR, confirming the asymptotic behavior observed in Fig. 9a 


3.3. Settling Fiber 

Consider a fiber settling under constant gravity force = Fj_ej_ acting perpendicularly 
to its major axis. The dynamics of the system depends on three competing effects: the 
elastic stresses which tend to return the object to its equilibrium shape, the gravitational 
acceleration which uniformly translates the object and the hydrodynamic interactions which 
creates local drag along the hlament. After a transient regime, the hlament reaches steady 
state and settles at a constant velocity with a hxed shape (see Figs. 10a and 10b). This 
steady state depends on the elasto-gravitational number 


B = F^L/K^ 


(65) 


173 HU and [65] examined the contribution of each competing effect by measuring the nor¬ 
mal deflection A, i.e. the distance between the uppermost and the lowermost point of the 
filament along the direction of the applied force (Fig. 


10 b) 


and the normal friction co¬ 
efficient 7 _l/ 75. as a function of B. 7 ° is the normal friction coefficient of a rigid rod. To 
compute hydrodynamic interactions m used Stokeslet ; jH], the Force Coupling Method 
(FCM) [IH] ; [65], Slender Body Theory. 


Similar simulations were carried out with both the Joint Model described in Section [2.6.2l 
and the Gears Model. Fiber of length L = 68 a is made out of = 31 beads with gap width 
2eg = 0.2a for the Joint Model and N}, = 34 for the Gears Model. To avoid both overlapping 
and numerical instabilities with the Joint Model, the following repulsive force coefficients 
were selected: do = (o + ^ 3 )/ 4 , Jd = (a + Ci = 0.01 and C 2 = 0.01. No adjustable 

parameters are required for Gears Model. 


Fig. 11 shows that our simulations agree remarkably well with previous results except 
slight differences with [63] in the linear regime B < 100. Using Slender Body Theory, [63] 
made the assumption of a spheroidal filament instead of a cylindrical one, with aspect ratio 
r„ = 100, i.e. three times larger than other simulations, whence such discrepancies. The 


normal friction coefficient (Fig. 11b ), resulting from hydrodynamic interactions, perfectly 
matches the value obtained by [H] with the Force Goupling Method. The FGM is known to 
better describe multibody hydrodynamic interactions. Such a result thus supports the use 
of the simple Rotne-Prager-Yamakawa tensor for this hydrodynamic system. 


Differences between Gears and Joint Models implemented here are quantified by mea¬ 
suring the relative discrepancies on the vertical deflection A 

Ag-Aj 


eyi = 


22 


(56) 












Figure 10: Shape of settling fiber for B = 10000 in the frame moving with the center of mass 
[xciZc). a) Metastable “W” shape, t = 12L/14- b) Steady “horseshoe” shape at t = 53L/Vs. 
Vg is the terminal settling velocity once steady state is reached. 


and on the normal friction coefficient 7 . 1 / 7 ° 


( 7 ^/ 7 o^)g- (7^/70^) j 

(7^/7o^)g 


(57) 


Discrepancies between Joint and Gears models remain below 5% except at the threshold 
of the non-linear regime {B ^ 100) where reaches 15% and ~ 7.5%. 


In accordance with ca. a metastable “W” shape is reached for B > 3000 (Fig. |10a[ ) 
until it converges to the stable “horseshoe” state (Fig. 10b). 


3.4- Actuated Filament 

The goal of the following sections is to show that the model we proposed is not only valid 
for passive objects but also for active ones. Elastohydrodynamics also concern swimming at 
low Reynolds number | 6 ]. Many type of micro-swimmers have been studied both from the 
experimental and theoretical point of view. Among them two categories are distinguished: 
ciliates and flagellates. Ciliates propel themselves by beating arrays of short hairs (cilia) 
on their surface in a synchronized way {Opalina, Volvox, Paramecia). Flagellates undulate 
and/or rotate their external appendage to push (pull) the fluid from their aft (fore) part 
(spermatozoa, Chlamydomonas Rheinardtii, Bacillus Subtilis, Eschericia Coli). Recent ad¬ 
vances in nanotechnologies allows researchers to design artihcial swimming micro-devices 
inspired by low Reynolds number fauna HlEgiTT]. 

In that scope, the study of bending wave propagation along passive elastic hlament has been 
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Figure 11: a) Scaled vertical deflection A/L depending on B. —•— ; Gears Model, () : 
Joint Model, -B— : FCM results from [H], -A— ; Stokeslets results from m,- -B— : Slender 
body theory results from [HS]- b) Normal friction coefficient vs. B. —•— : Gears Model, 
: Joint Model, -B— : FGM results from M.- -A- ; Stokeslets results from pg. 
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investigated by [TS] , [TU] and IHHl 121] ■ 


The experiment of [79] consists in a flexible filament tethered and actuated at its base. 
The base angle was oscillated sinusoidally in plane with an amplitude oq = 0.435rad and 
frequency 

Deformations along the tail result from the competing effects of bending and drag forces 
acting on it. A dimensionless quantity called the Sperm number compares the contribution 
of viscous stresses to elastic response [19] 


Sp 


K’’ j 


1/4 


L 


(58) 


7-*“ is the normal friction coefficient. When using Resistive Force Theory, is changed 

into a drag per unit length coefficient can be seen as the length scale at which bending 

occurs. Sp < 1 corresponds to a regime at which bending dominates over viscous friction: the 
whole filament oscillates rigidly in a reversible and symmetrical way. Sp^ 1 corresponds to 
a regime at which bending waves are immediately damped and the free end is motionless [T9] . 


The experiment of [HD] is similar to ra except that the actuation at the base is rota¬ 
tional. Here, the hlament was rotated at a frequency C forming a base angle Oo = 0.262rad 
with the rotation axis. 


In both contributions, the resulting hber deformations were measured and compared to 
Resistive Force Theory for several values of Sp. Simulations of such experiments |79|, 180] 
were performed with the Gears Model. 


3.4-1- Numerical setup and boundary conditions at the tethered base element 

Corresponding kinematic boundary conditions for BM are prescribed with the constraint 
formulation of the Euler-Lagrange formalism. 


Planar actuation. In the case of planar actuation di. we consider that the tethered, i.e. 
the hrst, bead is placed at the origin and has no degree of freedom 


' — 0 

Ujf = 0 . 


Denote ao the angle formed between e^, and 61^2- 
The trajectory of bead 2 must follow 

2a cos (tto sin ((t)) 
r^(t) = I 0 

2a sin (ao sin ((t)) 
25 


(59) 


(60) 




The translational velocity of the second bead f 2 (t) is thus constrained by taking the deriva- 

= 


tive of (60) 


—2aaoC cos {C,t) sin (oo sin (Cf)) 
0 

2aaoC cos (Ct) cos (cto sin (('f)) 


(61) 


Helical actuation. In the case of helical beating 
slightly off-centered with respect to the rotation axis 
inset). [21] measured a value (5o = 2mm with a hlament length varying from L = 2cm to 


, the anchor point of the hlament is 
r(0) = 5o (cf. Fig. Left 


10cm. Here we take (5o = 5osinao with Sq = 2.7a and vary the hlament length by changing 
the number of beads Nb to match the experimental range (5o/L = 0.1 —)■ 0.02. The position 
of bead 1 must then follow 


^0 cos (ao sin {(t)) cos (ao cos {(t)) 
= I ^0 cos (ao sin (Cf)) sin (tto cos (Ct)) 

^0 sin (ao sin (Cf)) 


(62) 


The translational velocity of the hrst bead fi (f) is thus constrained by taking the derivative 


of (62) 


m = 


( ^o«oC [“ cos (Cf) sin (ao sin (Cf)) cos (ao cos (Cf)) \ 

-|- sin (Cf) sin (ao cos (Cf)) cos (ao sin (C^))] 
(^ottoC [~ cos (Cf) sin (ao sin (Cf)) sin (ao cos(Cf)) 

— sin (Cf) cos (ao cos (Cf)) cos (ao sin (C^))] 
(5oaoC cos (Cf) cos (ao sin (Cf)) 

And the rotational velocity is set to zero (jJi = 0. 

The velocity of the second bead f 2 (^) is prescribed in synchrony with bead 1: 


V 


(63) 


/ 


^2{t) = 


/ (5o + 2a)aoC [~ cos (Cf) sin (ao sin (Cf)) cos (ao cos (Cf)) \ 
-|- sin (Cf) sin (ao cos (Cf)) cos (ao sin (C^))] 

(^0 + 2a)aoC [~ cos (Cf) sin (ao sin (Cf)) sin (ao cos (Cf)) 

— sin (Cf) cos (ao cos (Cf)) cos (ao sin (C^))] 

(^0 + 2a)aoC cos (Cf) cos (ao sin (Cf)) 


(64) 


The rotational velocity is consistently constrained by the no-slip condition. The three¬ 


dimensional curvature k is discretized with (34) 


In both cases, imposing actuation at the base of the hlament therefore requires the 


addition of three vectorial kinematic constraints, (59) and (61), to the no-slip conditions: 


A^c = 3(Aff, — 1) -|- 3 X 3. The additional Jacobian matrix writes 


ract 


I3 O3 O3 O3 
O3 I3 O3 O3 
O3 O3 I3 O3 
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(65) 












The corresponding right-hand side contains the imposed velocities 


-Qact _ 


0 

0 


( 66 ) 


for planar beating, and 


gac£ _ 


0 


(67) 


for helical beating. 

jact gact simply appended to ^ and B respectively; corresponding forces and 


torques are computed as explained before in Section 2.2 


3.4-2. Comparison with experiments and theory 

The dynamics of the system can be described by balancing elastic stresses (flexion and 
tension) with viscous drag. Subsequent coupled non-linear equations can be linearized with 
the approximation of small deflections or solved with an adaptive integration scheme (HU 


Planar actuation. |79] considered both linear and non-linear theories and included the effect 
of a sidewall by using the corrected RFT coefficients of [52] . 

Simulations are in good agreement with experiments, linear and non-linear theories for Sp 


= 1.73, 2.2, and 3.11 (Fig. 12). Even though sidewall effects were neglected here, the Gears 
Model provides a good description of non-linear dynamics of an actuated filament in Stokes 
flow. 

Helical actuation. Once steady state was reached, [50] measured the distance of the tip of the 


rotated filament to the rotation axis d = r{L) (cf. Fig. 13, Left inset). Figure 13 compares 
their measures with our numerical results. Insets show the evolution of the filament shape 
with Sp. The agreement is quite good. Numerical simulations slightly overestimate d for 
30 < Sp^ < 90. This may be due to the lack of information to reproduce experimental 
conditions and/or to measurement errors. As stated in [21], taking the anchoring distance 
^0 into account is important to match the low Sperm number configurations where Sq/L is 
non-negligible and the filament is stiff. If the anchoring point was aligned with the rotation 
axis (5o = 0), the distance to the axis of the rod free end would be d/L = sincto = 0.259 for 


small Sp, as shown on Fig. 13 


3.5. Planar swimming Nematode 

Locomotion of the nematode Caenorhabditis Elegans is addressed here as its dynamics 
and modeling are well documented [l3l|39|. C. Elegans swims by propagating a contraction 


wave along its body length, from the fore to the aft (Fig. 14a). modeling such an active 


hlament in the framework of BM just requires the addition of an oscillating driving torque 
7 '^(s,t) to mimic the internal muscular contractions. To do so, [H] used the preferred 
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Figure 12: Comparison with experiments and numerical results from na- Gears Model 
results are superimposed on the original Fig. 3 of | 71 ]. Snapshots are shown for four equally 

spaced intervals during the cycle for one tail with = 0.435rad. • : experiment, -- : 

linear theory, - : non-linear theory, - : Gears Model, a) C = 0.5 rad.s“^, Sp= 1.73. 

b) C = 1-31 rad.s“^, Sp= 2.2. c) C = 5.24 rad.s“^, Sp= 3.11. 
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with Sp^. Snapshots are shown for twenty eqnally spaced intervals during one period at 


steady state. Gray level fades as time progresses. Left inset: 5q/L is the distance of the 
tethered bead to the rotation axis, d/L is the distance of the free end to the rotation axis. 
(Main figure) Distance of the rod free end to the rotation axis normalized by the filament 
length d/L. () ■. experiment, • : Gears Model with no anchoring distance d^/L = 0, ■ : 
Gears Model with 5q/L = 0.1 —)■ 0.02 as in [21]. 
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curvature model. In this model, the driving torque results from a deviation in the centerline 
curvature from 

K^{s,t) =—Kq (s) sin {ks — 271 ft), (68) 

where Kq{s) is prescribed to reproduce higher curvature near the head: 




D 

0 


Ko, s < 0.5L 

2Kq (L — s) /L, s > 0.5L. 


(69) 


The amplitude Kq, wave number k and the associated Sperm number 


Sp = L (/p/A'‘) 


1/4 


(70) 


were tuned to reproduce the measured curvature wave of the free-swimming nematode. They 
obtained the following set of numerical values: Kq = 8.25/L, k = l.Svr/L and Sp* = 22.6^^^ 
. The quantity of interest to compare with experiments is the distance the nematode travels 
per stroke V/{fL). is assumed to be constant along s and is deduced from the other 


parameters. As for (32), the torque applied on bead i results from the difference in active 


bending moments across neighbouring links 


liif) = (si+i, t) - (si_i, t ), 


D 


D 


(71) 


with m'^ (si,t) = f)b(sj). is added to at step 3 of the algorithm in Section 

rai 

To match the aspect ratio of C. Elegans, Vp = 16, |13] put = 15 beads together 
separated by gaps of width 2eg = 0.2a. Here we assemble = 16 beads, avoiding the use 
of gaps, and employ the same target-curvature wave and numerical coefficient values. 


The net translational velocity V* = V/{fL) = 0.0662 obtained with our model matches 
remarkably well with the numerical results V/{fL) = 0.0664 |33] and experimental measure¬ 
ments V/{fL) Tn 0.07 [39] . 


3.6. Cooperative swimming 

One of the conhgurations explored in [H3] has been chosen as a test case for the interac¬ 
tions between in-phase or out-of phase swimmers. Two identical, coplanar C. Elegans swim 
in the same direction with a phase difference A0 which is introduced in the target curvature, 
and thus in the driving torque, of the second swimmer 

t) = —Kq ( s ) sin {ks — 27 t ft -|- A0). (72) 


The initial shape of the swimmers is taken from their steady state. We dehne d as the 
distance between their center of mass at initial time (see inset of Fig. 

Fig. 3, our results (Fig. 


14b). Similarly to 


14b) show that antiphase beating enhance the propulsion, 
whereas in-phase swimming slows the system as swimmers get closer. Even though the 
model swimmer here is different, the quantitative agreement with [HH] is strikingly good. 
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Figure 14: a) Simulated wave motion of a swimming model C. Elegans. The nematode 
swims leftward and gray level fades as time progresses. Motion is shown in a frame moving 
with the microswimmer center of mass, b) (Inset) Two C. Elegans beating in the same 
plane at a distance d in opposite phase (A0 = n). Nematodes swim leftward and gray level 
fades as time progresses. (Main hgure) Swimming speed of the center of mass of the system 
V normalized by the isolated swimming speed of C. Elegans V*. o ; in-phase motion 
(A0 = 0) ; □ : antiphase motion (A0 = tt). 
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Numerical work by [HI] also revealed that the average swimming speed of inhnite sheets in 
hnite Reynolds number flow is maximized when they beat in opposite phase. The conclusion 
that closer swimmers do not necessarely swim faster than individual ones has also been 
reported in [85]. They measured a decrease in the swimming speed of 25% for groups of 


house mouse sperm, as obtained on Fig. 14b for d/L = 0.2. 


3.7. Spiral swimming 

Many of the flagellate microorganisms such as spermatozoa, bacteria or artihcial micro¬ 
devices use spiral swimming to propel through viscous fluid. Propulsion with rotating rigid 
or flexible filaments has been thoroughly investigated in the past years ([HHl ESI EH ESI EH 
ISHj [77]). In this section we illustrate the versatility of the proposed model by investigating 
the effect of the Sperm number and the eccentricity of the swimming gait on the swimming 
speed of C. Elegans. 


3.7.1. Numerical configuration 

The target curvature of C. Elegans remains unchanged except that it is now directed 
along two components which are orthogonal to the helix axis. A phase difference A0 = 7r/2 
is introduced between these two components. The resulting driving moment writes: 

(sj, t) = aK^K^{si, t)ej_ -f l3K^K^{si, t, Afi = 7r/2)efe. (73) 


{e||,ex,ef,} are body fixed orthonomal vectors, ey is directed along the axis of the helix, 
is a perpendicular vector and is the binormal vector completing the basis (Fig. 15 


inset). The magnitude of the curvature wave along (resp. e?,) is weighted by a coefficient 
a (resp. fi). The trajectory of a body element in the plane {ex,efe} describes an ellipse 
whose eccentricity depends on the value of the ratio fi/a. When /3/q; = 0 the driving 
torque is two-dimensional and identical to the one used in Section 3.5 When fi/a = 1 the 


magnitude of the driving torque is equal in both direction, the swimming gait describes a 
circle in the plane {ex,ef,} (see Fig. 15 inset). For the sake of simplicity, here we take 

{e||,e_L,eb} = {e _ 

a = 1 and only fi is varied in the range [0; 1]. 


■x,ey,ez}. As in 3.4, the curvature is evaluated with (34). In the following. 


3.7.2. Results 


Figure lA compares the planar swimming speed of C. Elegans V*, with its “helical” 
version V, depending on the Sperm number defined in Section 3.5] (68) and on the ratio 
fi/a. The Sperm number Sp lies in the range [17^/'^; 1000^/^] = [2.03; 5.62]. The lower 
bound is dictated by the stability of the helical swimming. When Sp < 2.03, the imposed 
curvature reaches a value such that the swimmer experiences a change in shape which is not 
helical. This sudden change in shape breaks any periodical motion and makes irrevelant the 
measurement of a net translational motion. Such limitation is only linked to the choice of 
the numerical coefficients of the target curvature model. 

For the characteristic value Sp* = 22.6^/^ chosen by [43], the purely helical motion 
provides a swimming speed four times faster than planar beating. Even though the model 
swimmer is different here, this result qualitatively agrees with the observation of [56] for 
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Figure 15: Helical swimming of C. Elegans. (Inset) Snapshot for Sp* = 22.6^^^ and jS/a = 1. 
—•—; trajectory of the center of mass. (Main hgure) Swimming speed of the center of mass 

V normalized by the planar swimming speed of C. Elegans V*. - : (3/a = 1,- : 

(3/a = 0.5,-: 13/a = Q (planar motion), 

which spiral swimming was faster than planar beating. Beyond a critical value Sp ~ 2.6, 
planar beating is faster. For 13/a = 0.5 the swimming speed is always smaller than for planar 
beating except when Sp < 2.15. This last observation is not intuitive. A more extensive 
study on the effect of the eccentricity of the swimming gait on the swimming speed would 
be of interest. 

4. Conclusions 

We have provided a simple general theoretical framework for kinematic constraints to 
be used in three-dimensional BM. This framework permits to handle versatile and complex 
kinematic constraints between flexible assembly of spheres, and/or more complex non de¬ 
formable objects at low Reynolds numbers. Using Stokes linearity, this formulation requires, 
at each time step, the inversion of a 0{Nc x N^) linear system for an assembly having A/ 
constraints. Constraints are exactly matched (up to machine precision) and their evaluation 
is insensitive to time-step. Furthermore, since the formulation explicitly handles mobility 
matrices, it can be used with any approximation for hydrodynamics interactions, from free 
drain (no HI) to full Stokesian Dynamics. The proposed framework also implicitly incorpo¬ 
rates the generic influence of external flows on kinematic constraints, as opposed to previous 
BM formulation which necessitates some adjustments to the ambient flow in most cases. 

We also propose a simple Gears Model to describe flexible objects, and we show that such 
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model successfully predicts the fiber dynamics in an external flow, its response to an external 
mechanical forcing and the motion of internally driven swimmers. Qnantitative agreement 
with previons works is obtained for both slender objects (hbers, actnated hlaments) and 
non-slender swimmers {C. Elegans), allowing its nse in a wide variety of contexts. The 
Gears Model is easy to implement and fnlhlls several important improvements over previons 
BM : 


• There is no limitation on the hber cnrvatnre, since Gears Model does not need any 
repnlsive force nor gap width to be dehned. 

• Gears Model is more generic than previons ones, since there is no need for nnmerical 
parameter to be tnned. 

• When compared with Lagrange mnltiplier formnlation of Joint Model, Gears Model is 
also mnch more stable by two orders of magnitnde in time-step, a drastic improvement 
which offers nice prospects for the modeling of complex flexible assemblies. 

Finally it shonld be noted that even if we only consider simple collections of spheres, any 
complex assembly can be easily treated within a similar framework, which provide interesting 
prospects in the fntnre modeling of complex microorganims, membranes or cytoskeleton 
micro-mechanics. 
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Appendix A. Correspondence between M and ^ 


Matrix ^ dehned in (13) results from the rearrangement of the well-known mobility 


matrix M. This operation is necessary in order to combine constraint eqnation (|^ and 
mobility relation (12) to obtain the constraint forces ^c- 


Matrix M relates the collection of velocities V = (vi,..., v^r^) and rotations f2 = 
(lji, ..., to the collection of forces F = (fi,...,fArJ and torqnes F = ( 7 ;^,..., 7 ^J 


V 

n 


M^f yfiT 


F 

F 


(A.l) 


where is the 3Nb x 3Nb matrix relating all the bead velocities to the forces applied 
to their center of mass 




/ MYf 




Mjy. \ 


M'''' 

-LVl T 


YY / 
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(A.2) 






(A.l) is not consistent with the structure of the generalized velocity Q = (vi, lji, 
and force = (fi,7i, ...,f 7 V 6 , 7 ^J vectors. Thus we rearrange M into ^ such that 


ll 


M 


nr 


(A.3) 


to obtain a mobility equation suited for the Euler-Lagrange formalism 

/ Ql \ / ^ii ■ ■ ■ \ ^ 


(A.4) 


\ / V -^N^l ■ ■ ■ j V fA^i, / 


(A.4) is strictly equivalent to (13). 


Appendix B. Asymptotic limit of force and moment balance on the Gears 
Model 

In this appendix we show that slender body formulation for elastic fibers, when applied 
to Gears Model is consistent with the discrete formulation of force and moments balance 


(36) and (37) in the asymptotic limit of small beads. 

m 


The force balance equation for a beam is 

dn“ 


ds 


f = 0, 


(B.l) 


n“(s) is the resultant internal stress on a cross-section S{s) at arclength position s along 
the centerline 


n 


s = 


'S(s) 


t dS, 


(B.2) 


for which the tangent vector to neutral fiber centerline is also the unit normal vector to cross- 
section S{s). f is the force per unit length which contains any supplementary contribution 
to the internal elastic response of the material (e.g. hydrodynamic force per unit length). 
The moment balance reads 1551 


dm*' 

ds 


t X n“ + r = 0, 


(B.3) 


where m*^(s) is the moment of the flexion and torsion stresses on the cross-section which are 
related to the local deformation of Frenet-Serret coordinates and r is the torque per unit 
length resulting from supplementary contributions to the internal elastic response. 


Let us consider the curvilinear integral of (B.l) over each bead i, following the centerline 


of the skeleton joining the contact point Cj_i between bead i — 1 and i and the bead center 
Fj, as well as the bead center and the contact point Cj between bead i and z -|- 1 (see Figure 
1^. The curvilinear arclength s thus varies from 2ai to 2a{i -|- 1) within bead i. 
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Following the centerline, the integral of the internal stress contribution to (B.l) reads 


f2a{i+l) 

ds 


ds = n“(2a(z + l))-n“(2ai). 


(B.4) 


' 2ai 


In the limit of pointwise contacts, the normal stress produced by contact forces at contact 
point Cj_i located at point Xc._^ reads 




(B.5) 


where 6s stands for the surface Dirac distribution at the bead surface. Consequently the 
moment distribution associated with a Dirac contact forces applied at Xc^_^ is 




- Xc,_J = ati_i X 




(B.6) 


Since the area of the cross-section S{s) normal to the centerline tends to the bead surface 
itself as s —)■ 2ai or s —)■ 2a{i + 1), one can hnd that n*^(2ai) —?■ fc^_j and n®^(2a(z -|- 1)) —)■ 
as a —)■ 0 using (B.4), (B.5) and (B.2). Hence, the hnite size integral of (B.l), fulhlls the 

f = 0, (B.7) 


following limit as the bead radius tends to zero 


f - f 


which is consistent with the force used in (36). 

The second term of the moment balance equation (B.3) from contact point Cj_i to bead 
center r,- is 


-2a(i+l/2) 


t X n“ ds = tj_i X 


(‘2a{i+l/2) 


n“ ds = tj_i X 


' 2ai 


' 2ai 


'V,- 


o-dV] ■ ti_i. 


(B.8) 


where volume is the half-bead joining contact point Cj_i with bead center r^, whose 
pointing outward normal at r* is tj_i = The surface Si- enclosing half-bead Vi- is 

composed of half-sphere S^i- and disk ^j_. Si- = 5^i- U QJi-. Similarly, considering the 
moment balance equation (B.3) from bead center r* to contact point q leads to 

n2a{i+l) n2a{i-\-l) 


t X n“ ds = t,- X 


n“ ds = t,- X 


'2a(j+l/2) 


l2a{i+l/2) 


crdV ] ■ ti, 


'Vi+ 


(B.9) 


where volume is the half-bead joining bead center to contact point c,, whose pointing 
outward normal at is —t* = The surface Si+ enclosing half-bead Vi+ is composed 

of half-sphere <5^*+ and disk 5'^+ = U Hence the integrated contribution of 
the second term of the moment balance equation (B.3) is the sum of the right-hand-side of 
(B.8) and (B.9) which ought to be evaluated from the volume integral of the total stress over 
Vi- U inside bead i. Since the stress tensor inside the beads is not known, it is possible 
to relate it to the applied normal force at bead surface. Using divergence theorem on any 
volume V, enclosed by surface S, one finds 

/ aafsdV = / ■ n^)xfjdS = / cr • n 0 xc/F = D'^, 
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where n is the normal pointing outward surface S, whilst introducing the tensor as¬ 
sociated with the hrst moment contribution of the stress at surface S. If the surface S is 
the surface enclosing the considered bead, is the usual stress tensor, associated with 
the hydrodynamic interactions between the fluid and the bead. When considering hydrody¬ 
namic interactions, is usually decomposed into a symmetric tensor called stresslet and 


an anti-symmetric one called couplet. Using relation (B.IO) in (B.8) as well as (B.9), one 
hnds the following four contributions 

p2a{i-\-l) 

/ t X n“ ds = U _1 X + D^‘-) ■ ti_i + x (D-^*+ + D®*+) ■ t^, (B.ll) 

J 2a(i) 


to the integration of the second term of (B.3). In the limit of bead radius tending to zero. 


then tj_i —)■ tj, so that the outward normal vector to ^j_, tj_i, tends to the opposite of the 
outward normal vector to Since —>■ ^i+, this implies in turn, that —> —D^*+. 

Furthermore, since in the asymptotic limit of zero bead radius, U <5^,+ —)■ S, one hnds 
that 


/ t X n*® ds -)■ (tj_i X ■ ti_i) -F tj X ■ U)) -> U x (D"^ ■ U), (B.12) 

J 2a{i) 

where S is the bead surface here. For now, we concentrate on the contact forces contribution 
to (B.12). Using the contact surface force (B.5) it is then easy to compute the contact forces 
contribution to (B.12), 


r2a[i-\-l) 

/ t X n*^ds -)■ (ati_i x - aU x - rric,) 

J 2a{i) 


(B,13) 


Hence, result (B.13) is consistent with moment balance used in (37). 
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